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If magnetic frustration is most commonly known for undermining long-range order, as famously 
illustrated by spin liquids, the ability of matter to develop new collective mechanisms in order to 
fight frustration is no less fascinating, providing an avenue for the exploration and discovery of 
unconventional properties of matter. Here we study an ideal minimal model of such mechanisms 
which, incidentally, pertains to the perplexing quantum spin ice candidate Yb2Ti2 07. Specifically, 
we explain how thermal and quantum fluctuations, optimized by order-by-disorder selection, conspire 
to expand the stability region of an accidentally degenerate continuous symmetry U(l) manifold 
against the classical splayed ferromagnetic ground state that is displayed by the sister compound 
Yb2Sn207. The resulting competition gives rise to multiple phase transitions, in striking similitude 
with recent experiments on Yb2Ti2 07 [Lhotel et al., Phys. Rev. B 89 224419 (2014)]. Considering 
the effective Hamiltonian determined for Yb2Ti2 07, we provide, by combining a gamut of numerical 
techniques, compelling evidence that such multiphase competition is the long-sought missing key to 
understanding the intrinsic properties of this material. As a corollary, our work offers a pertinent 
illustration of the influence of chemical pressure in rare-earth pyrochlores. 


The wide interest in magnetic frustration largely 
stems from the diversity of unconventional phenomena 
it begets, ranging from anomalous Hall effect [l[ to mul- 
tiferroicity [H, to name only a few. The engine for this 
diversity is the inherent indecisiveness of frustrated mag¬ 
nets towards ordering, opening an avenue for new mech¬ 
anisms to control their low-temperature properties. The 
understanding of such mechanisms not only enlarges our 
broad understanding of the principles via which matter 
organizes itself; it is also the necessary first step in order 
to ultimately control the exotic properties of frustrated 
magnets. 


In the context of emergent degeneracies in py¬ 
rochlores ii, we present here a thorough analysis of 
the ordering mechanisms of a realistic microscopic model 
of multiphase competition. We believe our theory is the 
missing key to unlock many of the puzzling features of 
Yb2Ti207, whose ordering — or absence of — has been 


a matter of debate for nearly 15 yea^ 
the context of quantum spin liquid 
anism and magnetic monopoles 


, discussed in 
1 1211 . Higgs mech- 


14,113. Specifically, 


our theory accounts for the multi-step ordering process 
and field dependence of Yb2Ti207 [l 6 | while providing 
a natural setting to understand its sample dependence 


issue 


Starting from the most general model for anisotropic 
nearest- neig hbor exchange on the pyrochlore lat¬ 
tice [l 9 |, l 20 |, we consider a range of parameters known 
to describe Yb2Ti207 0 - We show that a splayed 
ferromagnetic (SFM) phase such as the one observed in 



FIG. 1. Multiphase competition of the anisotropic 

nearest-neighbor pyrochlore model of Eq. ©, as deter¬ 
mined by classical Monte Carlo simulations for {Jz=2,3,4} = 
{—0.22, —0.29, 0} meV. The -02 and phases are selected 
by order-by-disorder within the antiferromagnetic U(l) man¬ 
ifold and separated from the splayed ferromagnet (SFM) by 
a first-order transition, whose slope agrees with classical low- 
temperature expansion (green line, see Appendices). Rotating 
all spins of a given '02 state by the same angle around their lo¬ 
cal easy-axes produces the entire U(l) manifold, including the 
'03 states 0. The “splaying” angle of the spins in the SFM 
ground state is uniquely determined by the {Ji} parameters. 


the sister compound Yb2Sn2 07 [ 2 l|,[ 22 | competes against 
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a degenerate U(l) manifold for ordering (see Fig. [T]). As 
in Er 2 Ti 207 , a strong contender for experimental real¬ 
isation of order-by-disorder (ObD) |23l427|. soft modes 
of excitations lift the U(l) degeneracy. This degeneracy 
lifting creates an additional competition between the 
so-called ip 2 and states within the U(l) manifold. 
Using a wide range of techniques, we find that both 
thermal and quantum fluctuations enhance the stability 
of the U(l) manifold to the detriment of the SFM phase. 
As a consequence, quantum fluctuations bring Yb 2 Ti 2 07 
close to the boundary between all three phases where 
multiple phase transitions take place. In addition to the 
theoretical interest of understanding how multiple ObD 
selections take place within a classically energetically 
unstable manifold d, [HI, [HI , our work offers a natural 
explanation of the properties of Yb 2 Ti 207 , as recently 
observed in bulk measurements 0 - 


Model - We consider the anisotropic Hamiltonian 


n = Y,Siji 

M) 


'IJ 


with J = 


J2 Ja: 

— J 4 Ji J 3 
k T 4 1/3 J\ j 


( 1 ) 


where are pyrochlore nearest-neighbors. All coupling 
matrices Jij can be deduced from J by appropri¬ 
ate symmetry transformations 0 , with only four 
independent parameters {Ji=i.A} allowed by the sym¬ 
metry of the pyrochlore lattice 0 0 . Our work 
focuses on the parameter line Ji G [—0.09 : 0] meV 
and {Ji= 2 , 3 , 4 } = 1—0.22, —0.29, Q| meV, relevant to 
Yb 2 Ti 207 for Ji = —0.09(3) [s^ [i3 including 
the T = 0 boundary between SFM and U(l) phases 
at Ji = —0.029 (see Fig. [T]). We consider both quan¬ 
tum spins S = 1/2 and classical Heisenberg spins 
{\S\ = 1 / 2 ) whose phase diag rams have been studied 
in Refs. [ll|, [H, 26, 31, andIHj and Q, respectively. 


Classical thermal fluctuations - Starting at T = 0, the 
6 -fold degenerate SFM phase persists for Ji < —0.029 be¬ 
fore giving way to the one-dimensional degenerate U(l) 
manifold for Ji > —0.029. At the boundary, new contin¬ 
uously degenerate ground states emerge which confer ad¬ 
ditional zero-mode fluctuations to 7/3 configurations 0 . 

This zero-temperature framework sets the scene for 
the phase diagram of Fig. [1] computed by Monte Carlo 
(MC) simulations using parallel tempering [nisi and 
over-relaxation [HI- The U(l) degeneracy does not 
survive thermal fluctuations and collapses predomi¬ 
nantly in the 7/2 configurations, with a small 7/3 island 
around the boundary due to the above-mentioned soft 
modes of excitations. An interesting consequence of this 
ObD competition is that the ' 1 ^ 2/"^3 phases are precisely 
selected to optimize the entropy of the U(l) manifold 
for a given set of parameters and temperature. At finite 
temperature, this optimization puts the energetically 
selected SFM phase at a disadvantage and gives rise to 



Classical 

Quantum 


MC 

LSW 

ED 

NLC 

HTE 

T> 0 

-0.0340(5) 

nja 

nja 

-0.070(3) 

-0.06(3) 

T = 0 

-0.0289(1) 

-0.062 

-0.064(2) 

nja 

nja 


TABLE I. Critical value of the exchange parameter Ji [meV], 
separating the splayed ferromagnet (SFM) from the U(l) 
manifold, as estimated from different methods (Monte Carlo, 
linear spin-wave, exact diagonalization, numerical linked- 
cluster and high-temperature expansion) at zero temperature 
and upon cooling from high temperature. Quantum and ther¬ 
mal fluctuations work together to stabilize the U(l) manifold 
over the SFM phase. In particular, quantum effects bring the 
SFM/U(1) boundary within error bars of Yb 2 Ti 207 parame¬ 
ters pol] . 


multiple phase transitions for Ji G [—0.034 : —0.029]. 
Since the Hamiltonian of Eq. m supports a large variety 
of emergent degeneracies and potential ObD transitions 
at the boundaries between ordered phases [I, [HI, we 
expect such a phenomenology to be a common feature 
of pyrochlores [HI, 113 arid frustrated magnetism in 
general [HI- But since temperature is not the only 
source of fluctuations, how do quantum fluctuations fit 
in this picture ? 

Quantum fluctuations at zero temperature - Since we 
know the competing classical phases (SEM, ^2 and 7 / 3 ), it 
is of interest to analyze their stability in the semiclassical 
limit using linear spin-wave theory (LSW) [HI- However, 
when applied to classically unstable states, LSW usually 
becomes rather tedious as it requires higher order terms 
in the spin wave expansion. To circumvent this prob¬ 
lem, we used the method developed in Ref. [HI, which 
provides an upper bound of the semiclassical ' 02/'03 en¬ 
ergies for Ji < —0.029 meV (see Appendix). Keeping 
in mind that this approach underestimates the stability 
of the U(l) manifold, LSW shows that the semiclassical 
T = 0 frontier is shifted by quantum zero-point fluctu¬ 
ations from —0.029 meV down to —0.062 meV (see Ta¬ 
bled!). 

We now consider the full quantum nature of the spins. 
Since frustration is already at play in the constituting 
bricks of the pyrochlore lattice, namely the tetrahedra, 
exact diagonalization of a finite number of tetrahedra 
provides a good indication of the local influence of 
quantum effects. To preserve the symmetry of the 
pyrochlore lattice, we consider clusters of 4 and 16 spins, 
forming respectively 1 and 5 tetrahedra and allowing for 
standard ED calculations. Defining the order parameter 
M and correlator C = (M^) — (M)^ of a given phase, 
then the quantity AC = C'u(i) — Csfm is a direct 
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measure of the SFM/U(1) competition, bringing the 
T = 0 frontier to Ji —0.052(2) and —0.064(2) for 
N = 4 and N = 16 respectively, in agreement with the 
semi-classical results (see Table [Tj). 

Quantum fluctuations at finite temperature - Even if 
the joint analysis of thermal and quantum fluctuations 
is particularly challenging for such a frustrated prob¬ 
lem, the build up of correlations when approaching the 
phase transition from high temperature remains accessi¬ 
ble thanks to numerical linked-cluster computation 113 - 
lii and high-temperature expansion (see Appendix 
for technical details). 

HTE confirms the shifting of the boundary down to 
Ji = —0.06(3) meV and with transition temperatures 
lower than 500 mK (see inset of Eig. 0), in agreement 
with our classical simulations. As for NEC, at high 
temperature where quantum effects ultimately disap¬ 
pear, AC changes sign at the classical limit Ji ~ —0.03 
meV, which can be understood from I3‘^ terms in high- 
temperature expansion. Then, as temperature decreases 
for Ji < —0.03, instead of diverging towards SEM order¬ 
ing, AC shows a clear upturn towards enhanced U(l) 
correlations (see Eig. 0. This upturn is adiabatically 
evolving to lower temperature as Ji is decreased. NEC 
thus puts the U(1)/SEM frontier at Ji = —0.070(3) 
meV and, together with HTE, confirms the quantum 
nature of the boundary shift observed in ED and ESW, 
but now observed at non-zero temperature (see Table m. 


Multiphase competition in Yb 2 Ti^ 0^ - The Hamilto- 
nian o was shown to fit well inelastic neutron scatter¬ 
ing data of Yb 2 Ti 2 07 under large field 0 , allowing for 
the estimation of the exchange parameters {Ji=i, 2 , 3 , 4 } = 
{-0.09(3), -0.22(3), -0.29(2), 0.01(2)} meV. This set of 
parameters has been useful to understand the param¬ 
agnetic and high field regimes d, [l3, [H, [lH, El, 13 • 
But many questions remain open at low temperature 
and zero field. Eor instance, while single crystals show 
a broad peak at ~ 185 mK in specific heat, pow¬ 
der saimles typically display a sharp peak at ^ 265 
mK H, [ol, [3 5 suggestiM the influence of structural dis¬ 
order as in Tb 2 Ti 207 [3 or in Dy 2 Ti 207 |3- But in¬ 
terestingly, recent bulk measurements have brought to 
light ne w g eneric features between powder and single 
crystals [l6|. 

Based on these latter experiments [Til, thanks 
to the present analysis, we believe it is finally possible 
to flesh out a common framework for the powder and 
single crystals which magnetically order. In Ref. 0, 
both samples undergo a double phase transition with a 
high-temperature non-ferromagnetic transition followed 
by a low-temperature ferromagnetic first-order transi¬ 
tion. Eurthermore, the application of a magnetic field 
h does not destroy the lower-temperature ferromagnetic 
transition, but rather increases its temperature for h > 5 



FIG. 2. The difference of correlators AC — Gu(i) — Gsfm 
computed with NEC confirms the quantum shifting of the 
boundary towards more negative values of Ji than for the clas¬ 
sical system, estimated at Ji = —0.070(3) meV. Ji is given in 
the caption while {Ji= 2 , 3 , 4 } = {-0.22,-0.29,0}. Inset: AC 
as computed from HTE for Ji = —0.09 (red), —0.06 (cyan) 
and —0.03 (black). The breadth of each curve represents the 
uncertainty for NEC and HTE at low temperature. 


mT; the transition remaining first order up to he ~ 20 
mT, before becoming continuous or vanishing, which is 
experimentally difficult to distinguish. These experimen¬ 
tal results fit perfectly within the framework of our the¬ 
ory. The double transition is a direct consequence of the 
SEM/U(1) competition, as shown in Eig. [1] and for a simi¬ 
lar range of temperature as in experiments. Eurthermore, 
since the U(l) manifold is antiferromagnetic, it does not 
couple with h. The magnetic field thus only favors the 
SEM phase, and the first-order ferromagnetic transition 
is expected to persist until the U(l) phase is destroyed at 
he. Monte Carlo simulations confirm this scenario with 
he ~ 15 mT, followed by a crossover for h > he (see 
Fig. Hb). We also mention as a further similarity the 
experimental presence of a reversible bump in the sin¬ 
gle crystal magnetisation at ^ 180 mK, between the two 
transitions. According to our MC simulations, such fea¬ 
ture could correspond to a '02/'03 ObD transition, but 
with the caveat that the phase does not persist above 
50 mK in our classical phase diagram. In that case, strnc- 
tural disorder may play an important role 3 El, El , 
since it is known to i) favor over 7/^2 order [3 
ii) to be stronger in single crystals than powder samples 
where no bump was observed 0 . 


It should be noted that even if the parametrization 
that we used, taken from Ref. llOL was done on samples 
different from the ones in Ref. [l6|, the quantum shift¬ 
ing of the SEM/U(1) boundary would bring this classi¬ 
cal scenario within the experimental uncertainty of the 
Ji = —0.09 ± 0.03 meV parametrization range. This 
quantum shifting is best illustrated in Eig. O where the 
structure factor calculated from classical simulations at 
Ji = —0.0335 meV is almost identical to quantum (NEC) 
results at Ji = —0.06 meV. As for neutron scattering 
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FIG. 3. Application to Yh 2 Ti 2 07 : a) structure factor as measured by neutron scattering and b) phase diagram in a field, a) 
The spin-flip (SF, top panels) and non-spin-flip (NSF, bottom panels) calculated by quantum NLC and classical Monte Carlo are 
almost identical when Ji is shifted from —0.0335 meV to —0.06 meV (please note that the temperature has been renormalized 
by 5/4 for a better agreement). This agreement confirms the quantum shifting of the boundary at finite temperature. When 
approaching the phase transition, the MC structure factor for Ji = —0.0335 meV reproduces the characteristic features of 
Yb 2 Ti 2 07 neutron scattering data measured at 300 mK (to be compared with Fig. 2a of Ref. [l3|i. On the other hand, the 
comparison between experiments 0 and classical simulations using the parametrization of Ref. [lol| (Ji = —0.09 meV) are 
noticeably less successful, especially around (220) Q]. This disagreement is not a criticism of the parametrization by Ross et 
a/., but rather an emphasis on the importance of quantum fluctuations. The temperature of 0.65 K in the right panels has been 
chosen such that the ratio between measurement and transition temperatures is the same as in Ref. [l^ . The color scale is fixed 
from 0 to the maximum SF intensity, except for the panel at 0.3 K where the color scale was chosen for visual comparison with 
experiments 0. b) When a field is applied along the [001] direction, the antiferromagnetic ^2 phase gradually disappears in 
MC simulations. Our theory thus explains why the first-order transition only persists at low field in Yb 2 Ti 207 [3. We used 
Ji = —0.033 meV and the anisotropic g-tensor of Yb 2 Ti 2 07 0 : g± = 4.18 and pn — 1.77. 


measurements of Yb 2 Ti 207 (see Fig. 2.a of Ref. 0 ), 
the comparison is noticeably better with classical simu¬ 
lations for Ji = —0.0335 meV where a double transition 
takes place, than for J\ = —0.09 meV, confirming one 
more time the relevance of our theory to Yb 2 Ti 207 . 

Last but not least, our work brings Yb 2 Ti 207 as 
the missing link between Yb 2 Sn 2 07 and Yb 2 Ge 207 , 
whose ground states are respectively splayed fer¬ 
romagnetic [111, [111 and a yet not characterized 
antiferromagnet [49|, which we tentatively associate 
with U(l). On the basis of our work, we anticipate 
that a natural path will take this series of compounds 
through a transition from SFM to U(l), achieved either 
under chemical pressure (Sn-^Ti^Ge) or hydrostatic 
pressure, in analogy with spin ice [H^ and Tb 2 Ti 2 07 [Flj 
experiments, respectively. 

Conclusion - Using a palette of complementary nu¬ 
merical methods, our work sets a theoretical benchmark 
common to both powder and single crystals of Yb 2 Ti 207 . 
We have shown how the multi- step ordering recently ob¬ 
served in bulk measurements [l6| naturally arises from 
the competition between a SFM phase and a U(l) mani¬ 
fold mediated by order-by-disorder selection. On a more 
fundamental level, the underlying ordering mechanism 


should be understood as the joint consequence of both 
quantum and thermal fluctuations. The T = 0 frontier 
between the U(l) manifold and the splayed ferromagnet 
is shifted by quantum fluctuations, raising the interest¬ 
ing possibility to modulate frustration by tuning quan¬ 
tum fluctuations. In the general context of multiphase 
competition, order-by-disorder can be viewed as a free- 
energy optimization process which reinforces the stability 
of the degenerate phase it acts upon, by naturally se¬ 
lecting the subset of configurations with higher entropy 
and/or quantum zero-point fluctuations. 

In light of the numerous models and phases supported 
by pyrochlores, ranging from spin liquids and spin ices 
to (partially) ordered phases IllL l52l - 56j , and subsequent 
boundaries between them [^, [29|, [37 1, our present work 
is a paradigmatic example of why the properties of frus¬ 
trated magnets should generically be understood as the 
sum of competing phases, rather than coming from a sin¬ 
gle controlling state. Some of these properties would in¬ 
deed seem to be “coming from nowhere” in absence of a 
global phase diagram. We expect such competition be¬ 
tween neighboring phases to be particularly relevant to 
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proximity of phase boundaries. In that respect it is in¬ 
teresting to note that Er 2 Ti 207 , whose coupli ng p aram- 
eters lie far away from any phase boundary [sl, l2J, [261 , is 
one of the most robust rare-earth pyrochlore for repro¬ 
ducibility of experiments, while Yb 2 Ti 207 is essentially 
the antithesis. Hence, we expect the interplay between 
multiphase competition and disorder to become a very 
topical question, necessary to account for experiments in 
pyrochlores and frustrated magnetism, and promisingly 
rich in exotic physics. 

The authors would like to thank Han Yan for collabo¬ 
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search Chair program (M.G., Tier 1) and by the Perime¬ 
ter Institute (PI) for Theoretical Physics. Research at 
PI is supported by the Government of Canada through 
Industry Canada and by the Province of Ontario through 
the Ministry of Economic Development & Innovation. 

APPENDIX 

Order parameters and correlators 

As illustrated in Fig. 1 of the main text, three main 
phases are studied in this paper: the U(l) manifold, di- 

I 


vided into ^2 and phases, and the splayed ferromag- 
net (SFM). Their order parameters can be found in the 
literature d, [H, [H, ^ and are repeated below for 
convenience. 

In what follows, we adopt the convention of Ref. [l^ 
regarding the labeling of the spins, defined by their posi¬ 
tions relative to the centre of the tetrahedron 





The spins at the four sites are denoted as * 52 = 0 , 1 , 2,3 = 
{Sf, Sf, Sf) in the cubic coordinates, with \Si\ = 1/2. 

For the U(l) manifold, the order parameter, per tetra¬ 
hedron, is a two-dimensional vector. 


. \ 


mu(i) — 


m 


'U(l) 


V”^U(l)/ 


^ 1 f ^Qx \ cv \ cz qy Qz \ ^ Qx \ qy Qz \^qx qv \ Qz^ 

I ^^0 I *^0 I ^^2 I *^2 ^2 ' ^^3 *^3 ' ^ 3 ) 


V 


^^{-sy + s§^sy-s!-sy-si + sy + si) 


(3) 


/ 


To differentiate between the 1/2 and 1/3 states 
within the U(l) manifold, we define the angle 
^u(i) = arctan(my^^^/mg^^^); the function cos( 6 ^u(i)) 
respectively equals to -hi and —1 for 1/2 and 1/3 
states d, [ 2 §, [33- 


As for the splayed ferromagnetic phase, it is fully de¬ 
scribed by two three-dimensional order parameters d 


^SFMl = 


^SFMl 


/3 


^SFMl 


1 i 



+ Sf + + 5f) 

i{sy + + sy + Sf) 

+ Si + Si + SI) 


(4) 


^SFM2 = 


a \ 

^SFM2 

^SFM2 
\^^SFM2 j 


(5) 

\ 


—1 (qy I qz qy qz QV \ Qz _l_ qy qz\ 

^ Oq Di O2 -h 02 ^ O3 O3; 

—]-fqx _i_ qz qx i qz qx qz i qx qz\ 

-T Oq Di O2 02 -h O3 O3; 

— 1 / cx \ qy _ qx \ qy \ qx _ qy _ qx _ qy\ 

\2\/2'^*^0 ^ ^ -h 02 02 ^3 ^ 3 )J 


where msFMi is simply the global magnetization and 
hisFM 2 accounts for the splayed nature of the ferromag¬ 
netism, i.e. the fact that the spins are not colinear. For 
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a given SFM ground state, both msFMi and msFM 2 are 
finite. 

In order to compare the growth of correlations between 
the two phases, we define the order-parameter correlators 

Cij = {mimj) - {mi){mj) (6) 

where mi is a given order parameter. For the U(l) man¬ 
ifold, the use of the correlators is rather straightforward. 
However, for the SFM phase, since both msFMi and 
hisFM 2 are finite, one needs to compute the matrix of 
correlators 


C'SFM1,SFM1 C'sFMl,SFM2 
C'sFM2,SFMl C'sFM2,SFM2 


\ 


(7) 


By definition, this matrix is symmetric but a priori non¬ 
diagonal. Upon diagonalization, the maximum eigen¬ 
value is kept as correlator of the SFM phase. Please 
note that by symmetry of the lattice, calculations can 
be simplified by considering only one component of the 
order parameters, say e.g. ^spMi ^sfm 2 * 


Monte Carlo simulations 


The Fig. 1 of the main text has been obtained by 
classical Monte Carlo simulations of Heisenberg spins on 
the pyrochlore lattice. The conventional cubic unit cell 
consists of 16 spins and the system size has N = 3456 
spins. We used the standard Metropolis algorithm where 
a Monte Carlo step (MCs) is defined as N random single- 
spin-fiip attempts. To improve the quality of the simula¬ 
tions, parallel tempering [SSl, 341 and over-relaxation [35| 
were included in the simulations. 

The error bars of Fig. 1 (main text) were obtained by 
using two different cooling procedures during the equili¬ 
bration of the simulations. For the first “annealed” pro¬ 
cedure, the initial configuration is chosen randomly; the 
system is then gradually cooled down from high temper¬ 
ature (fixed at 10 K) to the temperature of measurement 
T during 10^ MCs; it is then equilibrated at the temper¬ 
ature T during 10^ additional Monte Carlo steps. Since 
it is starting from high temperature, the annealed pro¬ 
cedure tends to favor the phase with higher entropy, i.e. 
the U(l) manifold, and provides a lower boundary to 
the transition temperature. For the second “quenched” 
procedure, the initial configuration is fixed in the SFM 
phase; the system is then equilibrated at temperature T 
during 10^ additional Monte Carlo steps. Since it starts 
in the ordered SFM phase, the quenched procedure fa¬ 
vors the SFM phase and provides an upper boundary to 
the transition temperature. Following these equilibration 
procedures, measurements are done every 10 MCs during 
10^ MCs. 


For the phase diagram of Fig. 3 b) (main text), the 
system size was 3456 sites with measurements during 10^ 
MCs. For the structure factors of Fig. 3 a) (main text), 
no parallel tempering was necessary for simulations above 
the transition temperature. The system size was 128000 
sites with measurements during 10^ MCs. 


Classical low-temperature expansion 

Classical low-temperature expansion is an expansion 
in small fluctuations around an ordered state of classical 
spins. It enables calculation of the free energy of a given 
ordered phase up to leading term in temperature. We 
have used it to determine the low-temperature depen¬ 
dence of the phase boundary between the SFM and 
phases, for comparison with MC simulation as indicated 
by the green line in Fig. 1 of the main text. The method 
is a standard one, outlined in (e.g.) Ref. (H^ . 


Linear spin wave theory 

Linear spin wave theory (LSW) is a semi-classical 
method to study the stability of a classical phase in pres¬ 
ence of the quantum zero-point energy. When the phase 
is a classical ground state, such as the SFM phase in 
the double-transition region of our paper, the method is 
rather straightforward. But if the phase is not a classical 
ground state, such as the U(l) manifold in the double¬ 
transition region, the inclusion of zero-point energy re¬ 
quires a variation of the LSW theory, as proposed in 
Ref. [ 3 ^ and outlined below. 

At first, the approach is similar to standard spin wave 
expansion. We rewrite the spin operators Si in terms of 
Holstein-Primakoff bosons and perform a I/S' expansion 
around the local spin configuration of the chosen ordered 
state. At the harmonic ievei, our Hamiltonian takes the 

form n « "Hlsw = T^lsw + '^lsw + '^lsw- T^lsw is simply 
the classical energy of the ordered state around which we 

are expanding, while Hlsw ^lsw contain only linear 
and quadratic terms in boson operators, respectively. 

If the chosen ordered state is a classical ground state, 
then ^Lsw racist vanish. As for ^[1^, it can be diag¬ 
onalized by Fourier transform followed by a Bogoliubov 
transformation which will return real, positive, frequen¬ 
cies and therefore a meaningful excitation spectrum. The 
semi-classical energy can then be written 

= + ( 8 ) 

where is the classical ground state energy, is the 
spin wave dispersion of band A at wave-vector k. 

However, if the chosen ordered state is not a classi¬ 
cal ground state, then may still vanish if the spin 
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configuration we are expanding around is a saddle point 
of the classical energy. As shown in Ref. [ 2 ^, this is 

the case for the ^2 and configurations. The last re- 

( 2 ) 

maining point is the diagonalization of 'HWm' Following 
the approach developped by Coletta et al. [39| , we add to 
the positive definite term V = 5 ajoi, where the 


n 


( 2 ) 

LSW 


aj and Oi operators are Holstein-Primakoff bosons. This 
additional term V does not change the classical energy 
and the parameter 6 may be adjusted to the minimum 
value for which we can obtain a real, positive excitation 
spectrum. The energy calculated with the inclusion of V 




N- 


-Toji 

2 kX 

kX 


(9) 


where Ec\ is the classical energy of the saddle-point con¬ 
figuration and is the spin wave dispersion calculated 
including the potential V. Since R is a positive definite 
operator, F’semi-ci upper bound on the semiclassical 
energy of the saddle-point configuration. 


Exact diagonalization 

We consider the properties of the exchange Hamilto¬ 
nian shown in Eq. (1) of the main text for a single tetra¬ 
hedron {N = 4) and a single cubic unit cell of the py- 
rochlore lattice (A^ = 16) with periodic boundary con¬ 
ditions. Thanks to the small number of sites the full 
spectrum is analytically accessible for A" = 4 and the 


lowest lying eigenstates can be easily found via standard 
Lanczos diagonalization for A = 16. 


Numerical linked cluster expansion 

At finite temperature, we have also performed a nu¬ 
merical linked cluster expansion (NEC) 0, IIH, defined 
as 


P{C)/N = ^ L{C)W{C) (10) 

ccc 

where P is some extensive quantity and A is the number 
of lattice sites. The sum runs over clusters C of the lat¬ 
tice, L{C) counts the number of clusters of type C per site 
and W{C) is the weight evaluated on the cluster. This 
weight is computed using inclusion-exclusion rule 


W{C) = P{C) - Yi W{C') 

C'CC 


(11) 


where P(C) is the quantity computed on the cluster C and 
the sum runs over proper subclusters of C. There is some 
freedom in choosing the classes of clusters to sum over in 
this expansion. We follow the approach of Ref. and 
use tetrahedra as our building block. A linked cluster 
with ut tetrahedra will have at most Sut + 1 sites so we 
are limited to ut < 4 in the expansion. The properties 
P we will compute are defined on the tetrahedron and so 
conform well to this expansion. 


For convenience, we reproduce in Table [IT] the geomet¬ 
rical clusters used in Ref. , with the appropriate em¬ 
bedding constant. These include the 0^^ order point up to 
4 th which includes clusters composed of four tetra¬ 

hedra. Each cluster has an associated Hamiltonian He 
obtained from the exchange Hamiltonian (Eq. (1) of the 
main text) which we diagonalize numerically. The largest 
riT = 4 clusters we consider have 13 sites and thus Hilbert 
spaces of dimension 2^^ = 8192. These remain amenable 
to full diagonalization and thus we can compute arbitrary 
thermodynamic quantities at finite temperature. 


the embedding constants one has 


Po = PP{Co) 

Pi = -P(Co) + ip(Ci) 

P2 =-^PiCi) + P{C2) 

P 3 = +Ip{Ci) - 5 P(C 2 ) + 3P(C3) 

Pi = -Ip{Ci) + 10P(C2) - 21P(C3) 

+3P(C4a) + 6P(C45) -h 2P(C4 c) 


( 12 ) 

(13) 

(14) 

(15) 

(16) 


The series is organized into terms P^ including up to n 
tetrahedra. Explicitly carrying out the expansion using 


Note that P(Co) does not appear directly past first order. 
Following Ref. the Euler resummation method is used 
on the final two terms to accelerate convergence. If we 
define the differences Sn = Pn+i — Pn then the Euler 
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C L{C) 
Co 1 • 



C L{C) 


C2 



C4a 


C4c 



2 



TABLE IL Clusters used for the NLC expansion. A graphical representation is shown along with the embedding constant 

L{C). 


approximants are given by 

E 2 = P 2 (17) 

Es = E 2 ^ -Ss (18) 

^4 = ^3 + J (^3 + ^ 4 ) (19) 

The difference between 3^*^ and 4^^ order of expansion is 
the uncertainty of our NLC computations. 

High Temperature Expansion 

The High Temperature Expansions (HTE) are done up 
to order as explained in the book of Ref. [i^. The 
series expansion are then analyzed using Fade approxi¬ 
mants, i.e. based on rational functions. We constructed 
all near-diagonal Fade approximants with 8 or 7 terms 
in the series, i.e. [4/4], [5/3], [3/5], [6/2], [2/6], [4/3], 
[3/4], [5/2], [2/5] where [m/n] stand for the powers of 
the polynomial in the numerator and the denominator. 
The resulting spread in Fade approximant values repre¬ 
sents the error bars of our HTE computations. 
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